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ABSTRACT 

With the coming data dchigc from synoptic surveys, there is a growing need for frameworks that can 
quickly and automatically produce calibrated classification probabilities for newly-observed variables 
based on a small number of time-series measurements. In this paper, we introduce a methodology for 
variable-star classification, drawing from modern machine-learning techniques. We describe how to 
homogenize the information gleaned from light curves by selection and computation of real-numbered 
metrics (features), detail methods to robustly estimate periodic light-curve features, introduce tree- 
ensemble methods for accurate variable star classification, and show how to rigorously evaluate the 
classification results using c;ross validation. On a 25-class data set of 1542 well-studied variable stars, 
we achieve a 22.8% overall classification error using the random forest classifier; this represents a 
24% improvement over the best previous classifier on these data. This methodology is effective for 
identifying samples of specific science classes: for pulsational variables used in Milky Way tomography 
we obtain a discovery efficiency of 98.2% and for eclipsing systems we find an efficiency of 99.1%, both 
at 95% purity. We show that the random forest (RF) classifier is superior to other machine-learned 
methods in terms of accuracy, speed, and relative immunity to features with no useful class informa- 
tion; the RF classifier can also be used to estimate the importance of each feature in classification. 
Additionally, we present the first astronomical use of hierarchical classification methods to incorporate 
a known class taxonomy in the classifier, which further reduces the catastrophic error rate to 7.8%. 
Excluding low-amplitude sources, our overall error rate improves to 14%, with a catastrophic error 
rate of 3.5%. 

Subject headings: stars: variablcis: general - methods: data analysis - methods: statistical - tech- 
niques: photometric 
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1. INTRODUCTION 



Variable star science (e.g., E yer &: Mowlavi|2008 ) remains at the core of many of the central pursuits in astrophysics: 
pulsational sources probe stellar structure and stellar evolution theory, eruptive and episodic systems inform our 
understanding of accretion, stellar birth, and mass loss, and eclipsing systems constrain mass transfer, binary evolution, 
exoplanet demographics, and the mass-radius-temperature relation of stars. Some eclipsing systems and many of the 
most common pulsational systems (e.g., RR Lyrae, Cepheids, and Mira variables) are the fundamental means to 
determine precise distances to clusters, to relic streams of disrupted satellites around th e Milky Way, and to the local 
group of galaxies, 
recent review. 



They anchor the measurement of the size scale of the Universe. See Walkowicz et al. ( 2009 ) for a 



The promise of modern synoptic surveys ( Ivezic et al.||2007 ), such as the Large Synopti c Survey Telescope (LSST), 
is the promise of discovery of many new instances of variable stars (Sesar et al. 2007), some to be later studied 
individually with greater photometric and spectroscopic scrutinjj^ and some to be used as ensemble probes to larger 
volumes. New classes (with variability reflecting p hysics not previous ly seen) and rare instances of existing classes of 
variables are almost certainly on the horizon (e.g.. Covey et al.||2007 ). 

Classification of variable stars — the identification of a certain variable with a previously identified group ("class") 
of sources presumably of the same physical origin — presents several challenges. First, time-series data alone (i.e., 
without spectroscopy) provides an incomplete picture of a given source: this picture is even less clear the more poorly 
sampled the light curve is both in time and in precision. Second, on conceptual grounds, the observation of variability 
does not directly reveal the underlying physical mechanisms responsible for the variability. What the totality of the 
characteristics are that define the nature of the variability may in principle be known at the statistical level. But why 
that variability is manifest relies on an imperfect mapping of an inherently incomplete physical model to the data. (For 
example, the periodic dimming of a light curve may be captured with a small number of observable parameters but 
the inference that that source is an eclipsing one requires a theoretical framework.) This intermingling of observation 
and theory has given rise to a taxonomy of variable stars (for instance, defined in the GCV^ that is based on an 
admixture of phenomenology and physics. Last, on logistical grounds, the data volume of time-series surveys may be 
too large for human-intensive analysis, follow-up, and classification (which benefits from domain-specific knowledge 
and insight). 

While the data deluge problem suggests an obvious role for computers in classificatioij^ the other challenges also 
naturally lend themselves to algorithmic and computational solutions. Individual light curves can be automatically 
analyzed with a variety of statistical tools and the outcome of those analyses can be handled with machine-learning 
algorithms that work with existing taxonomies (however fuzzy the boundary between classes) to produce statistical 
statements about the source classification. Ultimately, with a finite amount of time-series data we wish to have 
well-calibrated probabilistic statements about the physical origin and phenomenological class of that source. 

While straightforward in principle, providing a machine-learned clas sifier that is accur ate, fast, and well-calibrated 
is an extraordinarily difficult task on many fronts (see discussion in Eyer et al. 2008 ) . There may be only a few 
instances of light curves in a given class ("labelled data") making training and vahdation difficult. Even with many 
labelled instances, in the face of noisy, sometimes spurious, and sparsely sampled data, there is a limit to the statistical 
inferences that can be gleaned from a single light curve. Some metrics (called, in machine-learning parlance, "features") 
on the light curve may be very sensitive to the signal-to-noise of th e data and others, particularly frequency-domain 
features, may be sensitive to the precise cadences of the survey (j |4.9[ ). For computationally intensive feature generation 
(e.g., period searches) fast algorithms may be preferred over slower but more robust algorithms 
Machine lear ning in variable star classification has been a pplied to several large time-scries datasets (Wozniak et al 



2004| jWillcmsen fc Eyer|[2007l [Debosscher et al.||2007] [Mahabal et al.n2008, ,Sarro et al.n2009^ |Blomme et al., 2010 



f 



A common thread tor most previous work is application of a certain machine-learning framework to a single survey. 
And, most often, the classification i s used to distin guish /identify a small set of classes of variables (e.g., Miras and 



other red giant variability). Debosscher et al. (^2007) was the first work to tackle the many-class (> 20) problem with 



multiple survey streams. [Debosscher et al. "(2007) also explored several classification frameworks and quantitatively 
compared the results. 

The purpose of this work is to build a many-class classification framework by exploring in detail each aspect of 
the classification of variable stars: proper feature creation and selection in the presence of noise and spurious data 
(Spj), fast and accurate classification (Sjsj), and improving classification by making use of the taxonomy. We present 
a formalism for evaluating the results of the classification in the context of expected statistical risk for classifying 
new data. We use data analyzed by Debosscher et al. to allow us to make direct comparison with those results (Q. 
Overall, we find a 24% improvement in the misclassification rate with the same data. The present work only makes 
use of metrics derivable from time-domain observations in a single bandpass; color information and context (i.e., the 
location of the variable in the Galaxy and with respect to other catalog sources) are not used. In future work, we will 
explore how machine-learned classifiers can be applied across surveys (with different characteristics) and how context 
and time-domain features can be used in tandem to improve overall classification. 



^ High-precision photometry missions (Kepler, MOST, CoRoT, etc.) are already challenging the theoretical understanding of the origin 
of variability and the connection of some specific sources to established classes of variables. 
^ General Catalog of Variable Stars, http://www.sai.insu.su/groups/cluster/gcvs/gcvs/ 
Not discussed herein are the challenges associated with discovery of variability. Sec Shin et al. i 2009 1 for a review. 
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2. HOMOGENIZING LIGHT CURVES: FEATURE GENERATION 

Classification fundamentally relies upon the ability to recognize and quantify the differences between light curves. 
To build a classifier, many instances of light curves are required for each class of interest. These labelled instances are 
used in the training and testing process (§31). Since the data are not, in general, sampled at regular intervals nor are 
all instances of a certain class observed witn the same number of epochs and signal-to-noise, identifying th e differences 
direct ly from the time-series data is exceedingly challenging both conceptually and computationally (cf. Eads et al 



2004 1 . Instead we homogenize the data by transforming each light curve into a set of real-number line features using 
statistical and model-specific fitting procedures. For variable stars, features fall into two broad categories: those 
related to the period of the source (and harmonics) and those that are not. Which features to use (and not use) is an 
important question that we will address herein. We also address the effects of (implicit) correlation of certain features 
in affecting the classification model. 

Appendix IX] provides an account of the non-periodic features used in this present work; many of these are simple 
statistics on the distribution of the fluxes (e.g., median absolute deviation and min-max amplitude) and some are 
domain specific (such as a feature that captures how much a source varies like the damped random walk seen in 
quasars; IButler & Bloom||2010l) . Since the period and periodic signatures of a variable are such crucial quantitative 
measurements, yet tend to be difficult to infer from simple prescriptions, we review the algorithms we employ to 
compute these features. 

2.1. Robust estimation of periodic features 
2.1.1. A Fast Period Search Including Measurement Uncertainty 

We model the photometric magnitudes of variable stars versus time t as a superposition of sines and cosines, starting 
from the most basic form: 

yz{t\fi) = a, sin(27r/jt) + hi cos{2tt fit) -f 6i,o, (1) 

where a and h are normalization constants for the sinuoids of frequency /j, and hi o is the magnitude offset. For each 
variable star, we record at each epoch, i^, a photometric magnitude, dk, and its uncertainty, Uk- To search for periodic 
variations in these data, we fit (fl]) by minimizing the sum of squares 



X 



Y^dk-y^(tk)Y/(^l 



(2) 



where is the measurement uncertainty in data point dk- As discussed in Zechmeister & Kiirster (20091, this least- 
squares fitting of sinusoids with a floating mean and over a range of test fre quencies is closely similar to an evaluation 
of the well-known Lomb-Scargle ( Lomb|[l976[ Earning 1963 Scargle 1982) periodogram. Allowing the mean to float 
leads to more robust period estimates in the cases where the periodic phase is not uniformly sampled; in these cases, 
the model light curve has a non-zero mean. (This is particularly important for searching for periods on timescales 
similar to the data span Ttot-) If we define: 



(3) 



with the weighted mean /i = ^k[dk/ '^l^t^ then our generalized Lomb-Scargle periodogram Pf is 



PfU) 



{N -i)xl- xlAf) 



xl 



(4) 



where Xmif) X^ minimized with respect to a, b, and 6o. For the NULL hypothesis of no periodic variation and a 
white noise spectrum, we expect Pf to be _F-distributcd with 2 nu merator an d — 1 denominator degrees of freedom. 
A similar periodogram statistic and NULL distribution is derived in Gregory) ( 2005 ) by marginalizing over an unknown 
scale error in the estimation of the uncertainties. In the lim it of many data, the NULL distribution takes the well- 
known exponential form (e.g., [Zechmeister fc Kurster||2009" I . For all = 1, Q becomes the standard Lomb-Scargle 
periodogram. In addition to the benefits of allowing a floating mean, the generalized Lomb-Scargle periodogram Q 
has two principal advantages over the standard formula: (1) uncertainties on the measurements are included, and 
because Pf{f) is a constructed from a ratio of two sums-of-squares, (2) there can be scale errors on the determination 
of these uncertainties (cf., ,Gregory„2005) . 

We undertake the search tor periodicity in each source by evaluating Q on a linear test grid in frequency from a 
minimum valu e of l/Ttpt to a maximu m value of 20 cycles per day, in steps of Sf = 0.1/Ttot- This follows closely the 
prescription in Debosscher et al. (|2007[), with the important exception that we search for periods up to 20 cycles per day 
in all sources, whereas 'Debosscher et al. (2007) search up to a "pseudo" Nyquist frequency (/jv — 0.5(1/Ar), where AT 
is the difference in time between observations and (•) is an average) for most sources but allow the maximum frequency 
value to increase for certain source classes. To avoid favoring spurious high-frequency peaks in the periodogram, we 
subtract a mild penalty of log ///jv above /at from P/(/) above / = /at. Significance of the highest peak is the evaluated 
from Pfif)- W e apply an approximate correction for the number of search trials using the prescription of Horne & 
Baliunas (1986), although we note that numerical simulations suggest these significance estimates underestimate the 
number oi true trials and are uniformly high by l-2cr. 




Fig. 1. — (A) Generalized Lomb-Scargle Periodogram Pf{f) for an eclipsing source in the sample. Plotted in blue is the first iteration 
to find peak frequency /i, which is twice the orbital frequency f-j (third iteration plotted in red). In this case, the second iteration yielded 
/2 = 3/3. For eclipsing sources, our Pf{f) analysis which utilizes a sine and cosine fit without harmonics, tends to place the orbital period 
in either /2 or fs. (B) The light curve folded at the orbital period /a. Over plotted is the best-fit model considering only /i (plus 3 
harmonics; in red), which fails to account for the difference in primary and secondary eclipse depths. Addition of the 2nd and 3rd frequency 
component models (black curve) account for the full light curve structure well. 



Press et al. 



Standard "fast" implementations of the Lomb-Scargle periodogram (e.g., 
number of frequency bins Nf as Nf log Nf, are not particularly fast for our purposes 
sample relatively few data (N ^ 100) on a very dense, logarithmic frequency grid Nf 

pursue algorithms which scale with N. We find that standard implementations are sped up by a factor of 



2001), which scale with the 
Tins is because we wish to 
10^. It is more fruitful to 
10 by 



simply taking care in calculating the sines and cosines that are necessary to tabulate Q. Instead of calculating all 
the sines and cosines at a given time point at each of the Nf frequency bins, we calculate sine and cosine only once at 
f = Sf and then use trigonometric identities (i.e., successive rotations by an angle 2^6 fti) to determine the sines and 
cosines at higher frequencies. 



2.1.2. Fitting Multiple Periods 



Following Debosscher et al. (2007), we fit each light curve with a linear term plus a harmonic sum of sinusoids: 



y(i) = rf + ^^2/,(t|i/,: 



(5) 



i=i j=i 



where each of 3 test frequencies fi are allowed to have 4 harmonics at frequencies fij = j fi. The 3 test frequencies fi 
are found iteratively, by successively finding the periodogram peaks in Pf(f). This procedure assumes a model with 
no harmonics, only the fund amental. We then r elax t his assumption by fitting for the fundamental plus 3 harmonics 
of the fundamental. 



Unlike Debosscher et al. 



(2007) who fit for the linear term prior to fitting sines and cosines, 
we evaluate the best-fit linear term at each test frequency during the search for the first periodogram peak, using a 
modification of (|4| to fit for the linear trend coefficient c at the same time as the constants ao, bo, and 60,0 • The value 
for c is updated when we then calculate the 4-harmonic fit (plus linear term) around the first best-fit frequency. The 



quantity Xo is updated after the i*^ iteration and subtraction of each 4-harmonic component is performed prior to 
calculating Pf{f)- 

In reporting the values from the fit of ([s]), we ignore the constant offsets 6^,0- We translate the sinusoid coefficients 
into an amplitude and a phase: 



A, 



Here Aj j (P Hj j) is the amplitude (phase) of the j**^ harmonic of the i"^ frequency component. Following 



(6) 
(7) 



Debosscher 



to relative phases with respect to the phase of the first component ^tl'i 



et al. (2007), we correct the phases PH^ j ^ ^..^^^ ^^^^^^^^ ^..^^ ...^ ^ — ^^..^..^ ^ ^^^ .^ 

Ptlij — PHqo ■ This is to preserve comparative utility in the phases for multiple sources by dropping a non-informative 
phase offset for each source. All phases are then remapped to the interval | — tt, -|-7r|. 
A list and summary of all of the period features used in our analysis is found in table [4] in Appendix 



2.2. Non-periodic light curve features 
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In seeking to classify variable star light curves, it may not always be possible to characterize flux variation purely 
by detecting and characterizing periodicity. We find that simple summary statistics of the flux measurements (e.g., 
standard deviation, skewness, etc.) — determined without sorting the data in time or period phase — give a great 
deal of predictive power. For instance, skewness is very effective for separating eclipsing from non-eclipsing sources. 
We define (in Appendix [A| 20 non-periodic features and explore in q4^thcir utility for source classification. A summary 
of all of the non-period features used in our analysis is found in tabic [5] in Appendix 

When only a small number of epochs (< 12) are sampled, we find that period detection becomes unreliable: we can 
only rely on crude summary statistics and contextual information to characterize these sources. In addition, some 
source classes yield non-period or multiply-periodic light curves, wherby the non-periodic features are expected in 
these cases to carry useful additional information not a lready contained in th e periodic features. As an example, we 
apply metrics used in the time-domain study of quasars ( Butler fc Bloom|2010 ) to aid in disentangling the light curves 
of some complexly varying, long-period variables (e.g., semi-regular pulsating variables). 



3. CLASSIFICATION FRAMEWORKS FOR VARIABLE STARS 

The features extracted from a light curve give a characterization of the observed astronomical source. We need a 
rigorous way of turning this information into a probabilistic statement about the science class of that source. This is 
the goal of classification: given a set of sources whose science class is known, learn a model that describes each source's 
class probabilities as a function of its features. This model is then used to automatically predict the class probabilities, 

and the most likely science class, of each new source. 

Seve ral authors have used machine-lear ning methods to class ify variable stars using their light curves: |Brett et al.| 
2004[) use Kohonen self-o rgani zing maps, Eyer fc Blake (2005) use the Bayesian mixture-model classifier Autoclass 



Cheeseman fc Stutz|1996 ) and Uebosscher et al. 12007) experiment with several methods, including Gaussian mixture 
models, Bayesian networks, Bayesian averaging ot artificial neural networks, and support vector machines. All of 
these methods have certainly enjoyed widespread use in the literature, and are a reasonable first set of tools to use in 
classifying variable stars. Our major contribution in this section is to introduce tree-based classification methods — 
including classification and regression trees (CART), random forests, and boosted trees — for the classification of 
variable stars. Tree-based classifiers are powerful because they are able to capture complicated interaction structure 
within the feature space, are robust to outliers, naturally handle multi-class problems, are immune to irrelevant 
features, easily cope with missing feature values, and are computationally efficient and scalable for large problems. 
Furthermore, they are simple to interpret and explain and ge nerally yield accurate res ults. In [|4]we show the superior 
performance of tree-based methods over the methods used in 'Debosscher et al. (2007) for classifying variable stars. 

Below, we describe several tree-based classification approaches from the statistics and machine learning literature, 
showing how to train each classifier and how to predict science-class probabilities for each observed source. We also 
introduce a suite of pair- wise voting classifiers, where the multi-class problem of variable-star classification is simplified 
into a set of two-class problems and the results are aggregated to estimate class probabilities. Additionally, we outline 
a procedure for incorporating the known variable star class taxonomy into our classifier. Finally, we describe a rigorous 
risk-based framework to choose the optimal tuning parameter(s) for each classifier, and show how to objectively assess 
the expected performance of each classifier through cross-validation. 



3.1. Tree-based classifiers 

Decision tree learning has been a popular m ethod for classification and regression in statistics and machine learning 



for more than 20 years (Breiman et al. 1984 popularized this approach). Recently, the astronomical community has 



begun to use t r ee-ba sed techniques, with impressive res ults. For exampl e, tre e-based classifiers hav e been used by 



Suchkov et al. 
separation, by 



( 20051) for SPSS o bject classification, by [Ball et al] ( |2006[ ) and [O'Keefe et al.| p009l ) for star-galaxy 
Bailey et ai.| (|2007p to identify supernova candidates, and by several grou ps torsupernova classification 



in the recent DFS Supernova Photometric Classification Challenge (Ke ssler et al.||2010 ) 



Tree-based learning algorithms use recursive binary partitioning to split the feature space, X , into disjoint regions, 
i?2, i?M- Within each region, the response is modeled as a constant. Every split is performed with respect 
to one feature, producing a partitioning of X into a set of disjoint rectangles (nodes in the tree). At each step, the 
algorithm selects both the feature and split point that produces the smallest impurity in the two resultant nodes. The 
splitting process is repeated, recursively, on all regions to build a tree with multiple levels. 

In this section we give an overview of three tree-based methods for classification: classification trees, random forest, 
and boosting. We focus on the basic concepts and a few particular challenges in using these cl assifiers for va riable-star 
classification. For further details about these methods, we refer the interested reader tO jHastie et al. (2009). 



3.1.1. Classification trees 

To build a classification tree, begin with a training set of (feature, class) pairs (Xi,yi), (X„,y„), where Yi can 
take any value in {1, C}. At node m of the tree, which represents a region i?,„ of the feature space X, the probability 
that a source with features in i?^ belongs to class c is estimated by 



Pr, 



- T 



(8) 
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which is the proportion of the Nm training set objects in node m whose science class is c, where I(Yi = c) is the 
indicator function defined to be 1 if 1^ = c and else. In the tree-building process, each subsequent split is chosen 
amongst all possible features and split points to minimize a measure of the resultant node impurity, such as the Gini 



index (X]^ c' PmcPmc' ) or entropy (— Pmc log2 Pmc) ■ The Gini index is the measure of choice for CART ( Breiman 



1996 ) . This splitting process is repeated 



observations in a terminal node, or 



et al.||1984[ ), while entropy is used by the popular algorithm C4.5 ( |Quinlan[ 

recursively until some pre-defined stopping criterion (such as minimum number o: 
relative improvement in the objective function) is reached. 

Once we have trained a classification tree on the examples (Xi, Yi), (X„, Y„), it is simple to ingest features from 
a new instance, X„ow, and predict its science class. Specifically, we first identify which tree partition X„cw resides in, 
and then assign it a class based on that node's estimated probabilities from ([s]). For example, if X„ow G i?„i, then the 
classification-tree probability that the source is in class c is 

Pc(Xnow)= Pmc (9) 

where Pmc is defined in (|8|. Using (9|, the predicted science class is p(Xnow) — arg nraXcPc(Xnow)- Note that we are 
free to describe the classmcation output for each new source either as a vector of class probabilities or as its predicted 
science class. 

There remains the question of how large of a tree should be grown. A very large tree will fit the training data well, 
but will not generalize well to new data. A very small tree will likely not be large enough to capture the complexity of 
the data- generating process. The appropriate size of a tree ultimately depends on the complexity of model necessary 
for the particular application at hand, and hence should be determined by the data. The standard approach to this 
problem is to build a large tree, T, with M terminal nodes and the n to prune this tree to find the subtree T* of T 



that minimizes a cross-validation estimate of its statistical risk (see ^3.5). 

3.1.2. Random forest 

Classification trees are simple, yet powerful, non-parametric classifiers. They work well even when the true model 
relating the feature space to the class labeling is complicated, and generally yield estimates with very small bias. 
However, tree models tend to have high variance. Small changes in the training set features can produce very different 
estimated tree structure. This is a by-product of the hierarchical nature of the tree model: small differences in the 
top few nodes of a tree can produce wildly different structure as those per turbations are p ropagated down the tree 



To reduce the variance of tree estimates, bagging (bootstrap aggregation, Breiman 1996 1 was proposed to average 
the predictions of B trees fitted to bootstrapped samples of the training data. Kanaom forest ( [Breiman_ 2001p is an 
improvement to bagging that attempts to de-correlate the B trees by selecting a random subset of mtry of the input 
features as candidates for splitting at each node during the tree-building process. The net result is that th e final 



averaged r andom forest model has lower variance than the bagging model, while maintaining a small bias (see Hastie 



et al.[ 2009, ch.l5 for a discussion). 



To obtain a random forest classification model, we grow B de-correlated classification trees. For a new variable star, 
the class probabilities are estimated as the proportion of the B trees that predict each class. As in classification trees, 
we are free to describe each source as a vector of class probabilities or a best-guess class. This prescription generally 
works well because by averaging the predictions over many bootstrapped trees, the estimated probabilities are more 
robust to chance variations in the original training set, and are almost always more accurate than the output of a 
single tree. Another advantage to random forest is the relative robustness of the estimates to choices of the tuning 
parameters (B, mtry; ^^nd the size of each tree) compared to other non-parametric classification techniques. In practice, 
we use the parameter values that give minimal cross-validation risk. 

3.1.3. Boosted trees 

Boosting is a method of aggregati ng simple rules to create a predictive model whose performance is 'boosted' over 
that of any of its ensemble members ( Freund & Schapire 1996) . In classification boosting, a sequence of simple classifiers 



(referred to as weak-learners) is applied, whereby in each iteration the training observations are re-weighted so that 
those sources which are repeatedly misclassified are given greater influence in the subsequent classifiers. Therefore, as 
the iterations proceed, the classifier pays more attention to data points that are difficult to classify, yielding improved 
overall performance over that of each weak-learner. The predicted class probabilities are obtained from a weighted 
estimate of the individual classifiers, with weights proportional to the accuracy of each classifier. 

Classification trees are natural base learners in a boosting algorithm because of their simplicity, interpretability, 
and ability to deal with data containing outliers and m issing values. M oreover, there arc efficient algorithms that can 
quickly estimate boosted trees using gradient boosting ( Friedman|2001 ). It is usually sufficient to use single-split trees 



(so-called decision stumps) as base-learners, though in situations with more complicated interactions, bigger trees are 
necessary. We use the training data to choose this tuning parameter through cross-validation. 

3.2. Measuring feature importance 

An additional advantage to tree-based classifiers is that, because the trees are constructed by splitting one feature at 
a time, they allow us to estimate the importance of each feature in the model. A feature's importance can be deduced 
by, for instance, counting how often that feature is split or looking at the resultant decrease in node impurity for splits 
on that feature. Additionally, random forest provides a measure of the predictive strength of each feature, referred 
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to as the variable importance, which tells us roughly what the decrease in overall classification accuracy would be if 
a feature were replaced by a random permutation of its values. Random forest has a rapid procedure for estimating 
variable importance via its out-of-bag samples for each tree, those data that were not included in the bootstrapped 
sample. 

Analyzing the feature importance is a critical step in building an accurate classification model. By determining 
which features are important for distinguishing certain classes, we gain valuable insight into the physical differences 
between particular science classes. Moreover, we can visualize which types of features are more predictive than others, 
which can inform the use of novel features or the elimination of useless features in a second-generation classifier. 



3.3. Pairwise classifiers 

A common approach for multi-class classification is to reduce the C-class problem into a set of C{C — l)/2 pairwise 
comparisons. This is a viable approach because two-class problems are usually easier to solve since the class bo undaries 
tend to be relatively simple. Moreover, some classification methods, such as support vector machines (SVM, Vapnik 
2000| , are designed to work only on two-class problems. In pairwise classification, classifiers for all C(C— 1)/2 two-class 
problems are constructed. The challenge, then, is to map the output from the set of pairwise classifiers for each source 
(pairwise probabilities or class indicators) to a vector of C class probabilities that accurately reflects the science class 
of that source. This problem is referred to as pairwise c oupling. 



1990 Friedman 



The simplest method of pairwise coupling is voting (Knerr et al. 
object is determined as the winner in a pairwise head-to-head vote, 
the pairwise class probabilities and tends to estimate inaccurate C-class probabilities 



1 



1996 |, where the class of each 

airwise voting is sub optimal because it ignores 
In situations where pairwise 



class probab ility estimates are available , voting is outperforme d by other methods such as the KuUback-Leibler based 



technique of Hastie & Tibshirani (1998) and the approaches of Wu et al. (20041, which reduce the problem to solving 



a linear system of equations. In this paper, we explore the use of both tree- based classifiers and SVMs i n pairwise 
classifi ers. To obtain C-class probabilities, we employ the second pairwise coupling method introduced by Wu et al^ 
(2004). We refer the interested reader to that paper for a more detailed description of the method and a review of 
similar techniques in the literature. 



3.4. Hierarchical classification 

In variable star classification, we have at our disposal a well-established hierarchical taxonomy of classes based on 
the physics and phenomenology of these stars and stellar systems. For instance, at the top level of our taxonomy, we 
can split the science classes into three main categories: pulsating, eruptive, and multi-star systems. From there, we 
can continue to divide the sub-classes until we are left with exactly one of the original 25 science classes in each node 
(see figure [2]). For classification purposes, the meaning of the hierarchy is clear: mistakes at the highest levels of the 
hierarchy are more costly than mistakes made at deeper levels because the top levels of the hierarchy divide physical 
classes that are considerably different, whereas deeper levels divide subclasses that are quite similar. 

Incorporating a known class hierarchy, such as that of figure [2] into a c lassification engine is a research field that 



has received much recent attention in the machine learning literature (see Silla & Freitas [2011 for a survey of these 
methods). By considering the class hierarchy, these classifiers generally outperform their fiat classifier' counterparts 
because they impose higher penalties on the more egregious classification errors. In th is paper, we consider two t ypes 
of hierarchical classification approaches: 
HMC (hierarchical multi-label classification, Blockeel et al 
forests of decision trees. Below, we provide a synopsis of 



2006) 



H SC (hierarchical sing le-label classification, Cesa-Bianchi et al 

y 



2006[ ) and 

We implement both HSC and HMC using random 
SC and HMC. For more details about these methods, see 



Vens et al. 



( |2008D . 

In HSC, a separate classifier is trained at each non-terminal node in the class hierarchy, whereby the probabilities 
of each classifier are combined using conditional pro babil ity rules to obtain each of the class probabilities. This has 
a similar flavor to the pairwise classifler approach of |3.3[ but by adhering to the class hierarchy it need only build a 
small set of classifiers and can generate class probabilities in a straight-forward, coherent manner. Moreover, different 
classifiers and/or sets of features can be used at each node in HSC, allowing for the use of more general cla ssifiers at 
the top of t he hierarchy and specialized domain-specific classifiers deeper in the hierarchy. A recent paper of |Blomme| 



et al. ( 2010| ) applied a method similar to HSC, using Gaussian mixture classifiers, to classify variable stars observed 
by the Kepler satellite. A second hierarchical classification approach is HMC, which builds a single classifier in which 
errors on the higher levels of the class hierarchy are penalized more heavily than errors deeper down the hierarchy. 
In the version of HMC that we use, the weight given to a classification error at depth d in the class hierarchy is Wq, 
where wq S (0, 1) is a tuning parameter. This force s the algorithm to pay more attention to the top level, minimizing 



the instances of catastrophic error (defined in ^4.3). 



3.5. Classifier assessment through cross-validation 

We have introduced a few methods that, given a sample of training data, estimate a classification model, p, to predict 
the science class of each new source. In this section, we introduce statistically-rigorous methodology for assessing each 
classifier and choosing the optimal classifier amongst a set of alternatives. Since our ultimate goal is to accurately 
classify newly-collected data, we will use the classifier that gives the best expected performance on new data. We 
achieve this by defining a classifier's statistical risk (i.e. prediction error) and computing an unbiased estimate of that 
risk via cross-validation. We will ultimately use the model that obtains the smallest risk estimate. 
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variable star 




Mira SRPV RVTau Class Popll MMode FM FO DM DScuti LEoo BCeph SPB GDor PVSC CP WR TTau HAEBE SDor Ell BPers BLyrae WUMaj 



Fig. 2. — Variable star classification hierarchy for the data used in !|4] This structure is a crucial element of the two hierarchical classifiers 
used in this study, HSC and HMC. The hierarchy is constructed based on knowledge of the physical processes that govern each type of 
object. At the top level, the sources split into three major categories: pulsating, eruptive, and multi-star systems. 

Given a new source of class Y, having features X, we define a loss function, L(Y,p(X)), describing the penalty 
incurred by the application of our classifier, p, on that source. The loss function encompasses our notion of how much 
the classifier p has erred in predicting the source's class from its features X. The expected value of L, E[L{Y,p{X.))], is 
the statistical risk, R{p), of the classifier. The expected value, E[-], averages over all possible realizations of (X,y) to 
tell us how much loss we can expect to incur for the predicted classification of each new source (under the assumption 
that new data are drawn from the same distribution as the training data). A key aspect to this approach is that 
it guards against over-fitting to the training data: if a model is overly complex, it will only add extra variability in 
classificatio n without decreasi ng the bias of the classifier, leading to an increase in the risk (this is the bias-variance 
tradeoff, see Wassermaii||2006 ) . 



Conveniently, within this framework each scientist is free to tailor the loss function to meet their own scientific 
goals. For instance, an astronomer interested in finding Mira variables could define a loss function that incurs higher 
values for mis-classified Miras. In this work, we use the vanilla 0-1 loss that is defined to be if F = p(X) and 1 if 
misclassified (here, p(X) is the estimated class of the source; alternatively, we could define a loss function over the set 
of estimated class probabilities, {pi(X), ...,pc(X.)}). Under 0-1 loss, the statistical risk of a classifier is its expected 
overall misclassification rate, which we aim to minimize. 

There remains the problem of how to estimate the statistical risk, R{p) of a classifier p. If labeled data were 
plentiful, we could randomly split the sources into training and validation sets, estimating p with the training set and 

computing a risk estimate R{p) with the validation set. Since our data are relatively small in number, we use fc-fold 
cross-validation to estimate R. In this procedure, we first split the data into K (relatively) equal-sized parts. For 

each subset k = 1,...,K the classifier is fitted on all of the data not in k and a risk estimate, i?(_fe)(p), is computed 
on the data in set k. The cross-validation risk estimate is defined as Rcvip) — 5= ^(-fc) (p)- shown in 



Burman 



cross- vahaation 



idati 



(1989), for K > 4, i?cv(p) is an approximately- unbiased estimate of R{p). In this paper, we use K — 10 fold 



In addition to selecting between different classification models, cross-validation risk estimates can be used to choose 
the appropriate tuning parameter(s) for each classifier. For instance, we use cross-validation to choose the optimal 
pruning depth for classification trees, to pick the optimal number and depth of trees to build, and number of candidate 
splitting features to use at each split for random forests, and to select the optimal size of the base learner, number 
of trees, and learning rate for boosted decision trees. The optimal set of tuning parameters for each method is found 
via a grid search. For each of the methods considered, we find that i?cv is stable in the neighborhood of the optimal 
tuning parameters, signifying that the classifiers are relatively robust to the specific choice of tuning parameters and 
that a grid search over those parameters is sufficient to obtain near-optimal results. 



4. CLASSIFIER PERFORMANCE ON OGLE+HIPPARCOS DATA SET 
4.1. Description of data 

In this paper, we test our feature extraction and classification methods using a mixture of varia ble star photometric 
data from the OGLE and Hipparcos surveys. OGLE (Optical Gravitational Lensing Experiment, Udalski et al.|[l999 l 
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TABLE 1 

Dataset characteristics, by survey. 



Survey 


NLC^ 


%NLCde6'^ 


NLC used<^ 


(rtot> (days)d (Ar,p„,h,> 


HIPPARCOS 


1044 


100.0 


1019 


1097 103 


OGLE 


523 


99.2 


523 


1067 329 



Total number of light curves available to us. 
^ Percentage of Debosscher et al. (20071 light curves available to us. 

Number of light curves atter removal of sources with ambiguous class and 
exclusion of small classes. 

Average time baseline. 



TABLE 2 

Dataset characteristics, by science class. 



Variable Star Class 


Namcdob^ 


NLC"" 


%NLCdei, 


Instrument 


{■^ epochs) 


min{fiY 






a. Mira 


MIRA 


144 


100.0 


HIPPARCOS 


98 


0.0020 


0.09 


11.2508 


b. Semireg PV 


SR 


42 


100.0 


HIPPARCOS 


99 


0.0010 


0.15 


1.0462 


c. RV Tauri 


RVTAU 


6 


46.2 


HIPPARCOS 


104 


0.0012 


0.05 


0.1711 


d. Classical Ccpheid 


CLCEP 


191 


97.9 


HIPPARCOS 


108 


0.0223 


0.15 


0.4954 


e. Pop. II Cepheid 


PTCEP 


23 


95.8 


HIPPARCOS 


107 


0.0037 


0.21 


0.7648 


f. Multi. Mode Cepheid 


DMCEP 


94 


98.9 


OGLE 


181 


0.5836 


1.21 


1.7756 


g. RR Lyrae, FM 


RRAB 


124 


96.1 


HIPPARCOS 


91 


1.2149 


1.95 


9.6197 


h. RR Lyrae, FO 


RRC 


25 


86.2 


HIPPARCOS 


92 


2.2289 


3.15 


4.3328 


i. RR Lyrae, DM 


RRD 


57 


100.0 


OGLE 


304 


2.0397 


2.61 


2.8177 


j. Delta Scuti 


DSCUT 


114 


82.0 


HIPPARCOS 


129 


0.0044 


7.90 


19.7417 


k. Lambda Bootis 


LBOO 


13 


100.0 


HIPPARCOS 


84 


7.0864 


12.36 


19.8979 


1. Beta Cephei 


BCEP 


39 


67.2 


HIPPARCOS 


96 


0.0014 


4.94 


10.8319 


m. Slowly Puis. B 


SPB 


29 


61.7 


HIPPARCOS 


101 


0.1392 


1.09 


11.8302 


n. Gamma Doradus 


GDOR 


28 


80.0 


HIPPARCOS 


95 


0.2239 


2.24 


9.7463 


o. Pulsating Be 


BE 


45 


78.9 


HIPPARCOS 


106 


0.0011 


2.12 


14.0196 


p. Per. Var. SG 


PVSG 


55 


72.4 


HIPPARCOS 


102 


0.0015 


3.41 


15.7919 


q. Chem. Peculiar 


CP 


51 


81.0 


HIPPARCOS 


105 


0.0076 


2.57 


13.4831 


r. Wolf-Rayet 


WR 


41 


65.1 


HIPPARCOS 


99 


0.0011 


6.56 


19.2920 


s. T Tauri 


TTAU 


14 


82.4 


HIPPARCOS 


67 


0.0013 


1.85 


11.2948 


t. Herbig AE/BE 


HAEBE 


15 


71.4 


HIPPARCOS 


83 


0.0009 


1.41 


10.0520 


u. S Doradus 


LBV 


7 


33.3 


HIPPARCOS 


95 


0.0008 


0.20 


0.5327 


V. Ellipsoidal 


ELL 


13 


81.2 


HIPPARCOS 


105 


0.1070 


1.37 


3.5003 


w. Beta Persei 


EA 


169 


100.0 


OGLE 


375 


0.0127 


0.93 


3.1006 


X. Beta Lyrae 


EB 


145 


98.6 


OGLE 


365 


0.0175 


0.71 


4.5895 


y. W Ursae Maj. 


EW 


58 


98.3 


OGLE 


369 


0.2232 


2.44 


8.3018 



^ Total number of light curves used, after removal of ambiguous sources. 
fl is the frequency of the first harmonic in day~^. 



is a ground-based survey from Las Ca mpanas Observatory cov ering fields in the Magellanic Clouds and Galactic bulge 



Hipparcos Space Astrometry Mission (Ferryman et al. [1997) was an ESA project designed to precisely measure the 
positions of more t han one hundred thousan d stars. 'Ifie data selected for this paper are the OGLE and Hipparcos 



sources analyzed by Debosscher et al. (2007), totaling 90% of the variable stars studied in that paper. A summary of 



the properties, by survey, of the data used in our study, is in table [T] The light curve data and classifications used for 
each source can he obtained thr ough our d ot astro . org light curve repository|^ 
This sample was designed by Debosscher et ah] ( |2007[ ) to provide a sizable set of stars within each science class. 



for a broad range of classes. Our sample contains stars from the 25 science classes analyzed in their paper. Via an 
extensive literature search, they obtained a set of confirmed stars of each variability class. In table [2] we list, by science 
class, the proportion of stars in that data set that we have available for this paper. Since the idea of their study was 
to capture and quantify the typical variability of each science class, the light curves were pre-selected to be of good 
quality and to have an adequate temporal sampling for accurate characterization of each science class. For example, 
the multi-mode Cepheid and double-mode RR Lyrae stars, which have more complicated periodic variability, were 
sampled from OGLE because of its higher sampling rate. 

In our sample, there are 25 objects that are labeled as two different science classes. Based on a literature search of 
these stars, we determine that 14 of them reasonably belong to just a single class (5 S Doradus, 2 Herbig AE/BE, 
3 Wolf Rayet, 2 Delta Scuti, 1 Mira, and 1 Lambda Bootis). The other 11 doubly-labeled stars, which are listed in 
table [6] were of an ambiguous class or truly belonged to two different classes, and were removed from the sample. See 
Appendix [B] fo r a detailed analysis and references for the doubly-labeled objects. Because the sample was originally 
constructeoby Debosscher et al. ( 2007 ) to consist only of well- understood stars with confident class labeling, we are 
justified in excluding these sources. 



' |http : //dotastro ■ org/llghtcurves/project ■php?Project_ID=123 
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TABLE 3 

Performance of Classifiers on OGLE+Hipparcos data set, 
averaged over 10 repetitions. 



Method Misclassification Catastrophic error %^ CPU*^ 



CART 


32.2 


13.7 


10.6 


C4.5 


29.8 


12.7 


14.4 


RF 


22.8 


8.0 


117.6 


Boost 


25.0 


9.9 


466.5 


CART.pw 


25.8 


8.7 


323.2 


RF.pw 


23.4 


8.0 


290.3 


Boost. pw 


24.1 


8.2 


301.5 


SVM.pw 


25.3 


8.4 


273.0 


HSC-RF 


23.5 


7.8 


230.9 


HMC-RF 


23.4 


8.2 


946.0 



^ Estimated using 10-fold cross-vaUdation. 

^ Average processing time in seconds for 10-fold cross-validation on 
a 2.67 GHz Macintosh with 4 GB of RAM. 



4.2. Class-wise distribution of light-curve features 

Using the methodology in ^ we estimate features for each variable star in the data set using their light curve. The 
feature-extraction routines talce 0.8 seconds per light curve, giving us a 53-dimensional representation of each variable 
star. The computations are performed in Python and C using a non-parallelized, single thread on a 2.67 GHz Intel 
Xeon X5550 CPU running on a v2.6.18 linux kernel machine. We estimate that the periodic-feature routines account 
for 75% of the computing time and scale linearly with the number of epochs in the light curve. Note that these metrics 
do not take into account the CPU time needed to read the XML data files from disk and load the data into memory. 

Plots of 1-dimensional density estimates, by science class, of a selected set of features are in figure|3] These class-wise 
feature distributions allow us to quickly and easily identify the differences, in feature space, between individual variable 
star science classes. Density plots are very useful for this visualization because they provide a complete feature- by- 
feature characterization of each class, showing any multi-modality, outliers, and skewness in the feature distributions. 
For instance, it is immediately obvious that several of the eruptive-type variable star classes have an apparent bi-modal 
or relatively flat frequency distributions, likely attributed to their episodic nature. Conversely, the RR Lyrae frequency 
distributions are all narrow and peaked, showing that indeed these stars are well-characterized by the frequency of 
their flux oscillations. The feature density plots also inform us of which feature(s) are important in separating different 
sets of classes. For example, the RR Lyrae, FO and RR Lyrae, DM stars have overlapping distributions for each of 
the features in figure |3] except the feature QSO, where their distributions are far apart, meaning that QSO will be a 
useful classification feature in separating stars of those two classes. 

4.3. Classifier comparison 

In this section, we compare the different classification methods introduced in SjSj To fit each classifler, except HMC- 
RF, we use the statistical environment if] To flt HMC-RF we use the open-source decision tree and rule learning 
system ClusFj Each of the classifiers was tuned via a grid search over the relevant tuning parameters to minimize 
the cross-valioation misclassification rates. To evaluate each classifier, we consider two separate metrics: the overall 
misclassification error rate and the catastrophic error rate. We define catastrophic errors to be any classification mistake 
in the top-level of the variable-star hierarchy in flgure[2](i.e., pulsating, eruptive, multi-star). The performance measures 
for each classifier, averaged over 10 cross-validation trials, are listed in table[3] In terms of overall misclassification rate, 
the best classifier is a random forest with B = 1000 trees, achieving a 22.8% average misclassification rate. In terms 
of catastrophic error rate, the HSC-RF classifier with B = 1000 trees achieves the lowest value, 7.8%. The processing 
time required to flt the ensemble methods is greater than the time needed to flt single tree models. However, this 
should not be viewed as a limiting factor: once any of these models is flt, predictions for new data can be produced 
very rapidly. For example, for a random forest classifler of 1000 trees, class probability estimates can be generated for 
new data at the rate of 3700 instances per second. 

4.3.1. Tree-based classifiers 

On average, the single-tree classiflers-CART and C4.5-are outperformed by the tree ensemble methods-random forest 
and tree boosting-by 21% in terms of misclas siflcation error rate. The classification performances of the single-tree 



classifiers are near that of the best classifier in Debosscher et al. (2007), which achieves a 30% error rate. Tree-based 
classifiers seem particularly adept at variable-star classification, and ensembles of trees achieve impressive results. 
Single trees take on the order of 10 seconds to fit a model, prune based on cross-validation complexity, and predict 
classification labels for test cases. Tree ensemble methods take 10-40 times longer to fit 1000 trees. Overall, the random 
forest classifler is the best classification method for these data: it achieves the lowest overall misclassification rate, low 
catastrophic error rate, and is the third fastest algorithm. 



^ R is a freely- available language and environment for statistical computing and graphics available at http://cran.r-project.org/ 
http://dtai.cs.kuleuven.be/clus/index.html 



Machine-Learned Classification of Variable Stars 



11 




QSO non_QSO 
Fig. 3. — Histograms of several features, by class. The features plotted are (a) the first frequency in cycles/day, (b) amplitude of the first 
frequency in mag, (c) the ratio of the second to the first frequencies, ( d) statistical significanc e of the periodic model , (e) flux ratio middle 
35th to middle 95th quantiles, (f) flux skew, (g) |Butler fc Bloom|201o| QSO variability, and (h) |Butler fc Bloom|2010| non-QSO variability. 
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4.3.2. Pairwise classifiers 

We implement four pairwise classifiers: CART, random forest, boosted trees, and SVM. Of these, random forest 
achieves the best results in terms of both misclassification rate and catastrophic error rate, at 23.4% and 8.0%, 
respectively. The pairwise classifiers all perform better than single-tree classifiers but tend to fare worse than the single 
random forest method. It is interesting to note that our implementation of SVM achieves a 25.3 % misclassification 
rate, a vast improvement over the 50% SVM misclassification rate found by Debosscher et al. (2007). This is likely due 
to both our use of better features and our pre-selection of the 25 features (chosen via cross-validation) with highest 
random-forest variable importance for use in the SVM. Unlike tree models, SVM is not immune to the inclusion of 
many useless features; when we include all 53 features into the SVM, our error rate skyrockets to 54%. 

4.3.3. Hierarchical classifiers 

Two of our classifiers, HSC-RF and HMC-RF, incorporate the hierarchical class taxonomy when building a classifier. 
Both of these methods achieve overall classification error rates slightly worse (sub-1% level) than that of the random 
forest classifier, while HSC-RF reaches the best overall catastrophic error rate (7.8%). HMC-RF slightly outperforms 
HSC-RF with respect to misclassification rate, but its current implementation takes 4 times as much CPU time. HSC- 
RF, HMC-RF, pairwise random forest and the original random forest are the best methods in terms of error rates, 
but random forest is at least twice as fast as any of the other methods. 



4.4. Direct comparison to Debosscher et al. {2007) 

Random forest achieves a 22.8% average misclassification rate, a 24% improvement over the 30% misclassification 
rate achieved by the best method of Debosscher et al. (2007) (Bayesian model averaging of artificial neural networks). 
Furthermore, each of the classifiers proposed in this paper, except the single-tree models CART and C4.5, achieves 
an average misclassification rate smaller than 25.8% (see table l3|. There is no large discrepancy between the different 
ensemble methods — the difference between the best (random forest) and worst (boosting) ensemble classifier is 2.2%, 
or an average of 34 more correct classifications — but in terms of both accuracy and speed, random forest is the clear 



In comparing our results to the Debosscher et al. (2007) classification results, it is useful to know whether the gains 
in accuracy are due to the use of better classifiers, more accurate periodic feature extraction, informative non-periodic 
features, or some combination of these. To this end, we classify these data using the following sets of features: 



1. The periodic features estimated by Debosscher et al. (2007). 



2. Our estimates of the periodic features following ^|2j 

3. The non-periodic features proposed in S|2] 

4. The non-periodic features in addition to our periodic features. 

Misclassification rates from the application of our classifiers to the above sets of features are plotted in figure |4] As 
a general trend, using both our periodic and non-periodic features is better than using only our periodic features, 
which is in turn better than using Debosscher et al.'s periodic features, which achieves similar rates to using only the 
non-periodic features. Using a random forest classifier, we find that the avera ge cross- validated misclas sification rates 



are 22.8% using all features, 23.8% using our periodic features, 26.7% using Debosscher et al. (|2007) features, and 



27.6% using only our non-periodic features. This is evidence that we obtain better classification results both because 
our classification model is better and because the extracted features we use are more informative. 



4.5. Analysis of classification errors 

Understanding the source and nature of the mistakes that our classifier makes can alert to possible limitations in the 
classification methods and feature extraction software and aid in the construction of a second-generation light-curve 
classifier. Let us focus on the classifier with best overall performance: the random forest, whose cross-validated mis- 
classification rate was 22.8%. In figure[5]we plot the confusion matrix of this classifier, which is a tabular representation 
of the classifier's predicted class versus the true class. A perfect classifier would place all of the data on the diagonal 
of the confusion matrix; any deviations from the diagonal inform us of the types of errors that the classifier makes. 

A few common errors are apparent in figure [5] For instance. Lambda Bootis and Beta Cephei are frequently 
classified as Delta Scuti. Wolf-Rayet, Herbig AE/BE and S Doradus stars are usually classified as Periodic Variable 
Super Giants, and T Tauri are often misclassified as Semiregular Pulsating Variables. None of these mistakes are 
particularly egregious: often the confused science classes are physically-similar. Also, the mis-classified exa mple s often 
come from classes with few training examples (see below) or characteristically low-amplitude objects (see ^4.8). 

As the random forest is a probabilistic classifier, for each source it supplies us with a probability estimate that 
the source belongs to each of the science classes. Until now, we have collapsed each vector of probabilities into an 
indicator of most probable class, but there is much information available to extract from the individual probabilities. 
For instance, in figure |6] we plot, by class, the random-forest estimated probability that each source is of its true 
class. We immediately see a large disparity in performance between the classes: for a few classes, we estimate high 
probabilities of true class, while for others we generally estimate low probabilities. This discrepancy is related to the 
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Fig. 4. — Distribution of cross-validation error rates for several classifiers on the OGLE-|-Hipparcos data, obtained from 10 repetitions. 
The classifiers are divided based on the features on which they were trained; from left to right: (1) all of the periodic and non-periodic 
features that we estimate, (2) the Lomb-Scargle features estimated by Debosscher ct al. (2007), (3) the Lomb-Scargle features estimated 
by us, and (4) the non-periodic features. In terms of mis-classification rate, the classifiers trained on all of our features perform the best, 
followed by those trained only on our periodic features, those trained on Debosscher et al.'s periodic features, and those trained on only 
the non-periodic features. All of the classifiers used, except single trees, achieve better error rates than the best classifier from Debosscher 
et al. (dashed line). 

size of each class: within the science classes that are data-rich, we tend to g et the correct clas s , whil e in classes with 



scarce data, we usually estimate the wrong class. This same effect is seen in Debosscher et al. (2007), (their table 5). 
This is a common problem in statistical classification for unbalanced class sizes: classihers such as random forest try to 
minimize the overall classification rate, thus focusing most of their efforts on the larger classes. One can remedy this 
by weighting the training data inversely proportional to their class size. Ho wever, the price to p ay for the increase in 



balanced error among the classes is a higher overall misclassification rate (see Breiman et al.|1984 ). The better solution 



is obviously to obtain more training examples for the under-sampled classes to achieve a better characterization of 
these classes. 

We have experimented with a random forest classifier using inverse class-size weighting. The results of this experiment 
were as expected: our overall misclassification rate climbs to 28.0%, a 23% increase in error over the standard random 
forest, but we perform better within the smaller classes. Notably, the number of correctly-classified Lambda Bootis 
increases from 1 to 7, while the number of correctly-classified Ellipsoidal variables jumps from 6 to 11, Beta Cephei 
from 5 to 23, and Gamma Doradus from 8 to 15. All four classes in which the original random forest found no correct 
classifications each had at least two correct classifications with the weighted random forest. 

4.6. Performance for specific science classes 

Although our classifier was constructed to minimize the number of overall misclassifications in the 25-class problem, 
we can also use it to obtain samples of objects from science classes of interest via probability thresholding. The random 
forest classifier produces a set of class-wise posterior probability estimates for each object. To construct samples of a 
particular science class, we define a threshold on the posterior probabilities, whereby any object with class probability 
estimate larger than the threshold is included in the sample. By decreasing the threshold, we trade off purity of the 
sample for completeness, thereby mapping out the receiver operating characteristic (ROC) curve of the classifier. 

In figure [Tj we plot the cross- validated ROC curve of the multi-class random forest for four different science classes: 
(a) RR Lyrac, FM, (b) T Tauri, (c) Milky Way Structure, which includes all Mira, RR Lyrae, and Cepheid stars, 
and (d) Eclipsing Systems, which includes all Beta Persei, Beta Lyrae, and W Ursae Major stars. Each ROC curve 
shows the trade-off between the efficiency and purity of the samples. At a 95% purity, the estimated efficiency for RR 
Lyrae, FM is 94.7%, for Milky Way Structure stars 98.2%, and for Eclipsing Systems 99.1%. The T Tauri ROC curve 
is substantially worse than these other classes due to the small number of sources (note: inverse class-size weighting 
does not help in this problem because the ordering of the posterior class probabilities drives the ROC curve, not the 
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Fig. 5. — Cross-validated confusion matrix obtained by the random forest classifier for the OGLE+Hipparcos data. A perfect classifier 
would place all mass on the diagonal. We have sorted the science classes by physical similarity, so large numbers near the diagonal signify 
that the classifier is performing well. On average, the classifier performs worse on the eruptive events, as exemplified by a large spread of 
mass for classes p through u. The overall error rate for this classifier was 21.1%. 

magnitude of those probabilities). Surprisingly, the 25-class random forest ROC curve dominates the ROC curve of a 
one-versus-all random forest for three of the four science classes, with vastly superior results for small classes. This 
implies that our approach of predicting all class probabilities with one classifier is better than constructing a separate 
classifier for each science class of interest. 



4.7. Feature importance 

In |3.2[ we described how to estimate the importance of each feature in a tree-based classifier. In figure [8] the 
importance of each feature from a pairwise random forest classifier is plotted, by class. The intensity of each pixel 
depicts the proportion of instances of each class that are correctly-classified by using each particular feature in lieu of a 
feature containing random noise. The intensities have been scaled to have the same mean across classes to mitigate the 
effects of unequal class sizes and are plotted on a square-root scale to decrease the influence of dominant features. In 
independently comparing the performance of each feature to that of noise, the method does not account for correlations 
amongst features. Consequentially, sets of features that measure similar properties-e.g. iiied.iaii_absolute_deviation, 
std, and stetson_j are all measures of the spread in the fluxes-may each have high importance, even though their 
combined importance is not much greater than each of their individual importances. 

Figure Is] shows that a majority of the features used in this work have substantial importance in the random forest 
classifier tor discerning at least one science class. Close inspection of this figure illuminates the usefulness of certain 
features for distinguishing specific science classes. As expected, the amplitude and frequency of the first harmonic, along 
with the features related to the spread and skew of the flux measurements, have the highest level of importance. The 
flux amplitude is particularly important for classifying Mira stars, Cepheids, and RR Lyrae, which are distinguished 
by their large amplitudes. The frequency of the first harmonic has high importance for most pulsating stars, likely 
because these classes have similar amplitudes but different frequencies. The QSO variability feature is important for 
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Fig. 6. — Cross-validated random forest probability estimate of the correct class for each variable star, plotted by class. Green circles 
indicate that the classifier's top choice for that object was correct, blue triangles indicate that the second choice was correct, and red x's 
indicated that neither of the top two were correct. As a general trend, we find that the classifier performs better for large classes. 

identifying eruptive sources, such as Periodically Variable Super Giants and Chemically Peculiar stars, because these 
stars generally have small values of the QSO feature compared to other variable star classes. 

In addition, there are several features that have very low importance for distinguishing any of the science classes. 
These features include the relative phase offsets for each of the harmonics, the amplitudes of the higher-order harmonics, 
and non-periodic features such as beyondlstd, max_slope, and pair_slope_trend. We rerun the random forest 
classifier excluding the 14 features with the smallest feature importance (the excluded features are 9 relative phase 
offsets, beyondlstd, max_slope, pair_slope_trend, and the 65th and 80th middle flux percentiles. Results show a 
cross-validated error rate of 22.8% and a catastrophic error rate of 8.2% (averaged over 10 repetitions of the random 
forest), which is consistent with the error rates of the random forest trained on all 53 features, showing the insensitivity 
of the random forest classifier to inclusion of features that carry little or no classification information. 

4.8. Classification of high- amplitude (>0.1 mag) sources 

High-amplitude variable stars constitute many of the central scientific impetuses of current and future surveys. In 
particular, finding and classifying pulsational variables with known period-lumi nosity relationships (Mira, Cepheids 



and RR Lyrae) is a majo r thrust of the Large Synoptic Survey Telescope (Walkowicz et al. 2009 'LSST Science 



Collaborations et al.||2009 ). Moreover, light curves from low-amplitude sources generally have a lower signai-to-noise 



ratio, making it more dihicult to estimate their light-curve derived features and hence more difficult to obtain accurate 
classifications. Indeed, several of the classes in which we frequently make errors are populated by low-amplitude 
sources. 

Here we classify only those light curves with amplitudes greater than 0.1 mag, removing low-amplitude classes from 
the sample. This results in the removal of 383 sources, or 25% of the data. After excluding these sources, we are left 
with a handful of classes with less than 7 sources. Due to the difficulty in training (and cross-validating) a classifier for 
classes with such small amount of data, we ignore those classes, resulting in the exclusion of 19 more sources, bringing 
the total to 408 excluded sources, or 26% of the entire data set. We are left with 1134 sources in 16 science classes. 

On this subset of data, our random forest classifier achieves a cross-validated misclassification rate of 13.7%, a 
substantial improvement from the 22.8% misclassification rate from the best classifier on the entire data set. The 




Fig. 7. — Cross-validated ROC curves, averaged over 10 cross validation trials, for four different science classes: (a) RR Lyrae, FM, (b) 
T Tauri, (c) Milky Way Structure, which includes all Mira, RR Lyrae, and Cepheid stars, and (d) Eclipsing Systems, which includes all 
Beta Persei, Beta Lyrae, and W Ursae Major stars. Plotted is 1-Efficiency versus 1-Purity as a function of the class threshold. The 25-class 
random forest ROC curve (solid black line) dominates the ROC curve of a one-versus-all random forest (dashed red line) for each science 
class except Eclipsing Systems. For T Tauri, the all-class random forest is vastly superior to the T Tauri-specific classifier. 

catastrophic misclassification rate is only 3.5%, compared to 8.0% for the entire data set. In figure [o] the confusion 
matrix for the classifier is plotted. The most prevalent error is misclassifying Pulsating Be and T Tauri stars as semi- 
regular pulsating variables. On average 8 of the 9 Pulsating Be stars and 4 of the 12 T Tauri stars are misclassified as 
semi-regular pulsating variables. 

4.9. OGLE versus Hipparcos 

Data from the sources that we classify in this section come from either the OGLE and Hipparcos surveys. Specifically, 
there are 523 sources from five science classes whose data are from OGLE, while the data from the remaining 858 
sources are from Hipparcos. Our classifiers tend to perform much better for the OGLE data than for the Hipparcos 
sources: for the random forest classifier we obtain a 11.3% error rate for OGLE data and 27.9% error rate for Hipparcos. 

It is unclear whether the better performance for OGLE data is due to the relative ease at classifying the five science 
classes with OGLE data or because of differences in the survey specifications. The sampling rate of the OGLE survey is 
three times higher than that of Hipparcos. The average number of observations per OGLE source is 329, compared to 
103 for Hipparcos, even though the average time baselines for the surveys are each near 1100 days. OGLE observations 
have on average twice the flux as Hipparcos observations, but the flux measurement errors of OGLE light curves tend 
to be higher, making their respective signal-to-noise ratios similar (OGLE flux measurements have on average a SNR 
1.25 times higher than Hipparcos SNRs). 

To test whether the observed gains in accuracy between OGLE and Hipparcos sources is due to differences in the 
surveys or differences in the science classes observed by each survey, we run the following experiment: For each OGLE 
light curve, we thin the flux measurements down to one-third of the original observations to mimic Hipparcos conditions 
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Fig. 8. — Pairwise random forest feature importance. Intensity is the square root of the proportion of instances of each class classified 
correctly because of that feature (compared to a replacement random permutation of that feature). Features are split into periodic (left) 
and non-periodic (right) features. The periodic features related to the first frequency are the most important, with higher-order frequencies 
and harmonics having smaller importance. The non-periodic features related to spread and skew have high importance in the classifier, as 
do the QSO-variability features. 

and re-run the feature extraction pipeline and classifier using the thinned light curves. Note that we do not add noise 
to the OGLE data because of the relative-similarity in average SNR between the surveys; the dominant difference 
between data of the two surveys is the sampling rate. Results of the experiment are that the error rate for OGLE 
data increases to 13.0%, an increase of only 1.7% representing 9 more misclassified OGLE sources. This value remains 
much lower than the Hipparcos error rate, showing that the better classifier performance for OGLE data is primarily 
driven by the ease of distinguishing those science classes. 



5. CONCLUSIONS AND FUTURE WORK 

We have presented a thorough study of automa ted variable star classifi cation from sparse and noisy single-band 
light curves. In the 25-class problem considered by 'Debosscher et al. (2007), which includes all of the most important 
variable star science classes, we obtain a 24% improvement over their best classifier in terms of misclassification error 
rate. Wc attribute this improvement to all of the following advances: 

• Better periodic feature estimation. Our Lomb-Scargle period-fitting code is both fast and accurate. With 
the same random forest classifier, the average error rate using our periodic feature estimates is 23.8%, compared 
to an error rate of 26.7% using only Debosscher's period feature estimates, representing an improvement of 11%. 

• Use of predictive non-periodic features. Simple summary statistics and more sophisticated model parame- 
ters give a significant improvement. Using both our periodic and non-periodic features, the random forest error 
rate is 22.8%, a 4% improvement over using only our periodic features. 
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Fig. 9. — Cross-validated confusion matrix for a random forest classifier applied only to the OGLE+Hipparcos sources with amplitude 
greater than 0.1 magnitude. The overall error rate for this subset of data is 13.7%, with catastrophic misclassification rate of 3.5%. 

• More accurate classification methods. All of the methods considered in this paper, save the single-tree 
models, achieve a statistically-significant improvement over Debosscher's best classifier. Our random forest 
classifier, applied to the exact features used by that paper, achieves an 11% improvement over their best error 
rate, 30%. 

Another contribution is our discussion of model validation through estimation of the expected prediction error on 
new data by cross-validation. We presented cross-validation as a statistically-rigorous way to both tune a classifier 
and select between competing methods. Indeed, all of the numbers quoted in this paper are 10-fold cross-validation 
estimates. 

We have shown the adeptness of tree-based classifiers in the problem of variable star classification. We demonstrated 
the superiority of the random forest classifier in terms of error rates, speed, and immunity to features with little useful 
classification information. We outlined how to calculate the optimal probability threshold to obtain pure and complete 
samples of specified sub-classes, and showed that the multi-class random forest is often superior to the one-versus-all 
random forest in this problem. We advocate the continued use of this method for other classification problems in 
astronomy. 

Furthermore, we described how the random forest classifier can be used to estimate the importance of each feature 
by computing the expected classification gains versus replacing that feature with random noise. In the variable star 
classification problem, it was found that several non-periodic features have high importance. A classifier built only on 
the non-periodic features still performs quite well, attaining 27.6% error rate. 

Finally, this paper is the first to use the known variable-star taxonomy both to train a classifier and evaluate its 
result. We introduced two different classification methods to incorporate a hierarchical taxonomy: HSC, which builds 
a different classifier in each non-terminal node of the taxonomy, and HMC, which fits a single classifier, penalizing 
errors at smaller depths in the taxonomy more heavily. We demonstrated that both of these methods perform well, in 
terms of classification rate and catastrophic error rate. The class taxonomy was also used to construct the notion of 
catastrophic error rate, which considers as catastrophic any error made at the top level of the hierarchy. 

Several open questions remain with regard to the automated classification of astronomical time series. Many of these 
questions will be addressed by us in future publications, where we will expand on the methodology presented here and 
attempt to classify data from other surveys, such as the All Sky Automated Survey (AS AS), SDSS Stripe 82, and the 
Wide Angle Search for Planets (WASP). Some of the questions that we will address are: 

• If we train a classifier on a set of objects from one (or multiple) survey (s), will that classifier be appropriate to 
predict the classes of objects from the new survey? This question is of great importance because presumably a 
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set of known (labeled) variable stars will be compiled from previous surveys to train a classifier for use on a new 
survey. 

• What features are robust across a wide range of different surveys, each with different cadences? If some sets of 
features are robust to survey design and cadence, those should be used in lieu of survey-dependent features in 
a classifier. In this paper, we have excluded any feature that was blatantly survey-dependent (such as any that 
used mean flux), but this does not guarantee that some features will not have survey dependence. 

• How does mis-labelled training data affect the classifier accuracy? Can mis-labelled data be effectively detected 
and cleaned from the classification set? 

• How can a classifier be trained to efficiently identify outliers/new types of variables? Future surveys will unques- 
tionably discover new science classes that do not fit under any of the training classes. Recent methodology has 
been developed in the statistics literature for outlier discovery in random forest classifiers, and we plan to adapt 
this for variable star classification. 

• How are the error rates of a classifier affected by computational limitations (where, perhaps some CPU-intensive 
or external server-dependent features are not used)? In automated classification of astronomical sources, there 
is often a time sensitivity for follow-up observations. Presumably there are more useful features for classification 
than the ones that we employed in this paper, but they may be expensive to compute or retrieve for each 
observation. This trade-off between error rate and computation time must be explored. 

Finally, as a longer-term goal, we are striving to develop methodology that can be used on a LSST-caliber survey. 
This means that our methods must be fast enough to compute features and class probabilities for thousands of objects 
per night, work well at an LSST cadence, be applicable to multi-band light curves, and perform classification for all 
types of astronomical objects, including transients, variable stars, and QSOs. Our task, looking forward, is to address 
each of these problems and develop methodology for fast and accurate classification for LSST. 
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APPENDIX 

A. FEATURES ESTIMATED FROM LIGHT CURVES 

A description of the 32 periodic features computed using the methodology in j]2.1|is in table [4l In addition to these 
32 periodic features, we calculate 20 non-periodic features for every light curve (j j2^ . These features are compiled 
in table [sj These consist primarily of simple statistics that can be calculated in the limit of few data points and also 
when no period is known in order to characterize the flux variation distribution. Where possible, we give the name 
of the Python function that calculates the feature (e.g., skewness is from scipy . stats . skewO in the Python SciPy 
module) . 

We begin with basic moment calculations using the observed photometric magnitude mag vector for each source: 

• skew: skewness of the magnitudes: scipy.stats.skew(); 

• small_kurtosis: small sample kurtosis of the magnitudej^ 

• std: standard deviation of the magnitudes: Numpy std(); 

• beyondlstd: the fraction (< 1) of photometric magnitudes that lie above or below one std() from the weighted 
(by photometric errors) mean; 



• stetson_j: Stetson (19961 variability index, a robust standard deviation; 

• stetson k: 



Stetson ( 1996 1 robust kurtosis measure. 



We also calculate the following basic quantities using the magnitudes: 

• max_slope: examining successive (time-sorted) magnitudes, the maximal first difference (value of delta magni- 
tude over delta time); 

• smiplitude: difference between the maximum and minimum magnitudes; 



see http: //wwM.xycoon. com/peakedness_small_sample_test_l .htm 
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median_absolute_deviation: mediandmag — median(TOa(7)|); 



• median_buf f er_r£m.ge_percentage: fraction (< 1) of photometric points within ampHtude/10 of the median 
magnitude 

• pair_slope_trend: considering the last 30 (time-sorted) measurements of source magnitude, the fraction of 
increasing first differences minus the fraction of decreasing first differences. 



We also characterize the sorted flux = 10 distribution using percentiles, following Eyer (2006). If ^5,95 is 

the difference between 95% and 5% flux values, we calculate: 



flux. 


.percentile. 


.ratio. 


.mid20: 


ratio ^40,6o/-p5,95; 


flux. 


.percentile. 


.ratio. 


.mid35: 


ratio F32.5,67.5/-F5,95 


flux. 


.percentile. 


.ratio. 


.midSO: 


ratio -F25,75/-P'5,95; 


flux. 


.percentile. 


.ratio. 


.mid65: 


ratio Fi7.532.5/-F5,95 


flux. 


.percentile. 


.ratio. 


.midSO: 


ratio -Fio,9o/-p5,95; 



• percent_amplitude :the largest absolute departure from the median flux, divided by the median flux; 

• percent_diff erence_f lux_percentile : ratio of -F5.95 over the median flux. 



Fina lly, useful for stochastically varying sources, we calculate the quasar similarity metrics from |Butler fc Bloom] 
(IMol): 



• QSQ: quality of fit Xqso/^ ^'^'^ ^ quasar-like source, assuming mag = 19; 

• non_QSO: quality of fit for a non-quasar-like source (related to NULL value of QSO). 



B. DOUBLY-LABELED STARS 



As discussed in |4.1[ there are 25 sources with more than one class label. Five of these objects are labeled as both S 
Doradus and periodically- variable super giants. We assign these to the S Doradus class because S Doradus is a subclass 
of periodicall y-variable super gia nts. Two other sources, V * BF Ori and HD 97048 were verified to be Herbig AE/BE 
type stars by Grinin et al. (20101 and Doering et al. (2007), respectively. RY Lep is listed incorr ectly in Simbad as an 
Agol-type eclipsmg system (its original 50-year old classitication) and as RR Lyrae/Delta Scuti. Derekas et al. ( 2009| ) 
clearly establi sh this source as a high-ampltiude J-Scuti in a binary system. Likewise, HIP 99252 is a d-bcuti star 
(Rodriguez et al. 2000) and not also a n RR Lyrae. HIP 30326 (=V Mon) is a Mira variable and not also part of the 
RV 'i'auri class (Whitelock et al.||2008 ). HD 210111 (listed also as delta Scuti) is a A Bootis star ( Paunzen & Reegen 
I. HD 1657BT is a WR star that is not periodic to a high-degre e of confidence (iMoffat et al.||2008p. Likewise, 
is not periodic as it was once believed ( Marchenko et aL]|1994 ). We found no reference to the periodic nature 
of the WR star HD 156385. The remaining 11 sources were of ambiguous class or truly deserved of two classes, and 
were excluded from the training and testing sample. These stars are listed in table |6] and briefly mentioned below: 



2008 
WW. 



• HD 136488 = HIP 75377: This is a WR with a measured periodicity in Hipparcos data ( |Koen fc Eyer||2002 1 

• HD 94910 ^ AG Car. This is a WR is thought to be a "hot, quiescent state LBV" ( [Clark et al^|2010[ ). 



• HD 50896 = EZ CMa = WR 6: Listed as both a WR and periodic va riable super giant, this is a prototype WR 
has a known periodic variability possibly due to a binary companion (Flores et al. 112007). 



HD 37151: Slowly pulsating B-stars (SPBs) occur in the same instability strip of the H-R diagram as chemically 
peculiar B stars (with varia bility associated wit h rotation and magnetic fields) . The classification of this star is 
particularly ambiguous (see Briquet et al.|[2007 and references therein). 



HD 35715 Is a cataloged Beta Cepheid but is also an ellipsoidal variable ( Stankov fc Handler||2005 ) . Pulsation 
was not detected photometrically but in line profiles. 

HIP34O42 (Z CMa A) Is a. binary system consisting of a Herbig component (A) and an FU Orionis star (B) in a 
close orbit. Determining the photometric contributions of the two components separat ely is problematic because 



of strong variability and the small angular separation ( Millan-Gabet fc Monnier 



parat e. 
[2002t . 



TVLyn, UY Cam, V1719 Cyg (HD 200925) & V753 Cen Ambiguous classification (within the literature) of RR 
Lyrae or Delta Scuti type stars. 
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TABLE 4 

Periodic features extracted from light curves using generalized Lomb-Scargle 



Feature 


Description'' 


f reql_harmonics_amplitude_0 




f reql_harmonics_amplitude_l 


Al,2 


f reql_harinoiiics_amplitude_2 


^1,3 


f reql_harinoiiics_amplitude_3 


Al,4 


f reql_harmoiiics_f req_0 


/l" 


f reql_harmoiiics_rel_phase_0 




f reql_harmoiiics_rel_phase_l 




f reql_harmoiiics_rel_phase_2 


PHl,3 


f reql_harmoiiics_rel_phase_3 


PHi,4 


f req2_harmoiiics_amplitude_0 


M,! 


f req2_harinoiiics_amplitude_l 


A2,2 


f req2_harmoiiics_amplitude_2 


^2,3 


f req2_harinoiiics_amplitude_3 


^2,4 


f req2_harinoiiics_f req_0 


h 


f req2_harmoiiics_rel_phase_0 


PH2,l 


f req2_harinoiiics_rel_phase_l 


PH2,2 


f req2_harmoiiics_rel_phase_2 


PH2,3 


f req2_harmoiiics_rel_phase_3 


PH2A 


f req3_harinoiiics_amplitude_0 


A3,l 


f req3_harinoiiics_amplitude_l 


^3,2 


f req3_harmoiiics_amplitude_2 


^3,3 


f req3_harinoiiics_amplitude_3 


^3,4 


f req3_harmoiiics_f req_0 


/3 


f req3_harmoiiics_rel_phase_0 


PH3,l 


f req3_harmoiiics_rel_phase_l 


PH3,2 


f req3_harmoiiics_rel_phase_2 


PH3,3 


f req3_harmoiiics_rel_phase_3 


PH3A 



f req_sigiiif Significance of /i vs. null hypothesis of white noise with no periodic 

variation, computed using a Student's-T distribution 
f req_signif _ratio_21 Ratio of significance of /2 vs. null to /i vs. null 

f req_sigiiif _ratio_31 Ratio of significance of /s vs. null to /i vs. null 

f req_varrat Ratio of the variance after, to the variance before subtraction 

of the fit with /i and its 4 harmonics 
f req_y_of f set c 

^ Notation from discussion of Lomb-Scargle periodic feature extraction in §2.1| is used. 
^ All amplitudes are in units of magnitude. 

All frequencies are in units of cycles/day. 

All relative phases are unitless ratios. 



TABLE 5 

Non-periodic features extracted from light curves 



Feature 



Description 



amplitude 
beyondlstd 

f lux_percentile_ratio_mid20 
f lux_percentile_ratio_mid35 
f lux_percentile_ratio_mid50 
f lux_percentile_ratio_mid65 
f lux_percentile_ratio_mid80 
linear_trend 
mcLx_slope 

mediaii_absolute_deviation 
median_buf f er_range_perceiitage 
pair _slope_t rend 
per cent .amplitude 

per cent _difference_flux_per cent ile 
QSO 

non_qSO 
skew 

smal 1 _kur t o s i s 
std 

stetson_j 
stetson_k 



Half the difference between the maximum and the minimum magnitude 
Percentage of points beyond one st. dev. from the weighted mean 
Ratio of flux percentiles (60th - 40th) over (95th - 5th) 
Ratio of flux percentiles (67.5th - 32.5th) over (95th - 5th) 
Ratio of flux percentiles (75th - 25th) over (95th - 5th) 
Ratio of flux percentiles (82.5th - 17.5th) over (95th - 5th) 
Ratio of flux percentiles (90th - 10th) over (95th - 5th) 
Slope of a linear flt to the light curve fiuxes 

Maximum absolute fiux slope between two consecutive observations 

Median discrepancy of the fluxes from the median flux 

Percentage of fiuxes within 20% of the amplitude from the median 

Percentage of all pairs of consecutive fiux measurements that have positive slope 

Largest percentage difference between either the max or min magnitude and the median 

Diff. between the 2nd & 98th flux percentiles, conv erted to magnitude^ 

Quasar variability metric in Butler fc Bloorn] ||20l"o[l 

Non-quasar variability metric m |Liutler 6c J:ilobm| ( |:2010^ 

Skew of the fiuxes 

Kurtosis of the fluxes, reliable down to a small number of epochs 
Standard deviation of the fluxes 
Welch-Stetson variability index j'' 
Welch-Stetson variability index K'' 



Iyer] |[2005ll 



Stetson ( 1996 i 
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TABLE 6 

Doubly-labelled sources excluded from the sample for both training and testing purposes. 



RA 



Dec 



ID 



Class 1 



Class 2 



RF best class 



RF 2nd best class ' 



05 26 

05 36 

06 54 

07 03 
07 33 
07 58 

10 56 

11 51 

12 04 
15 24 
21 04 



50.2284 
06.2333 
13.0441 
43.1619 
31.7288 
58.8801 
11.5763 
15.3088 
47.2738 
11.3087 
32.9211 



+03 05 
-07 23 
-23 55 
-11 33 
+47 48 
+72 47 
-60 27 
-55 48 
-27 40 
-62 40 
+50 47 



44.428 
47.320 
42.011 
06.209 
09.823 
15.411 
12.815 
15.795 
43.295 
37.567 
03.276 



HD 35715 
HD 37151 
HD 50896 
HIP 34042 
TV Lyn 
UY Cam 
HD 94910 
V753 Cen 
IK Hya 
HD 136488 
V1719 Cyg 



1. Beta Ccphei 
q. Chem. Peculiar 
r. Wolf-Rayet 
t. Herbig AE/BE 
j. Delta Scuti 
j. Delta Scuti 
r. Wolf-Rayet 
j. Delta Scuti 
e. Pop. II Cepheid 
p. Per. Var. SG 
j. Delta Scuti 



Ellipsoidal 
Slowly Puis. B 
Per. Var. SG 
T Tauri 
RR Lyrae, EG 
RR Lyrae, EG 
S Doradus 
RR Lyrae, EG 

g. RR Lyrae, EM 
r. Wolf-Rayet 

h. RR Lyrae, EG 



m. Slowly Puis. B (0.214) v, 

p. Per. Var. SG (0.211) n, 

V. Ellipsoidal (0.181) o. 

b. Semireg PV (0.287) d, 

h. RR Lyrae, EG (0.686) g. 

h. RR Lyrae, EG (0.236) j. 

b. Semireg PV (0.258) a. 

h. RR Lyrae, EG (0.469) j. 

g. RR Lyrae, EM (0.456) y. 

r. Wolf-Rayet (0.161) p, 

j. Delta Scuti (0.315) y. 



Ellipsoidal (0.164) 
Gamma Doradus (0.155) 
Pulsating Be (0.145) 
Classical Cepheid (0.114) 
RR Lyrae, EM (0.12) 
Delta Scuti (0.222) 
Mira (0.243) 
Delta Scuti (0.151) 
W Ursae Maj. (0.100) 
Per. Var. SG (0.147) 
W Ursae Maj. (0.226) 



Classes determined by [Debosscher et aT] | |2007| l . 

Most likely class determined by a random forest classifier. Posterior probability in parentheses. 
Second-best candidate class determined by a random forest classifier. Posterior probability in parentheses. 

• IK Hya Is a short-period Pop. II Cepheid, making it hkely a BL Hercuhs type star. 

We trained a random forest classifier on the features of only the single-labeled stars, and then predicted the class 
probabilities of each label for the doubly-labeled stars. This experiment shows that for 7 of the 11 sources, our classifier 
predicts one of the two classes that the star was originally labeled as. Note that recently, several method s have been 
developed to classify data with multiple labels (multi-label classification, see Tsoumakas fc Katakis||2007 1 , but for the 
purposes of this work we choose to only perform single-label classification, and remove these data from our sample to 
avoid biases from inaccurate training labels. 
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